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ABSTRACT 

Recent  research  by  Harrison  on  the  stability  of  parallel  flows  re- 
sulted in  a  successful  solution  of  plane  Poiseuille  flow  but  produced 
unexplained  anomalies  for  pipe  flow.   The  purpose  of  the  research  in 
this  paper  was  to  find  and  correct  errors  in  Harrison's  initial  analysis 
of  the  pipe  flow  problem. 

A  minor  error  in  Harrison's  numerical  method  was  corrected.  More- 
over, it  proved  necessary  to  reanalyze  the  conditions  on  the  axis  of 
symmetry  of  the  pipe.   These  changes  finally  made  it  possible  to  obtain 
reasonable  results. 

Owing  to  time  limitations,  the  number  of  solutions  obtained  using 
the  corrected  program  was  sufficient  only  to  confirm  its  general  vali- 
dity. However,  the  results  obtained  are  significant  in  that  they  dis- 
close instabilities  which  are  known  to  exist  but  which  have  not  been 
accounted  for  in  previous  theoretical  investigations. 
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2 

D,D  ,...   The  partial  derivatives  with  respect  to  r  in  cylindrical 

coordinates. 


e        Base  of  natural  logarithms.  Exponentiation  is  also  de- 
noted by  exp(  ) . 

e  ,  e  ,eQ     Unit  vectors  along  the  x,r,  and  Q  axes  in  cylindrical 

X.   J?   w 

coordinates. 

F,G,H     Components  of  the  velocity  vector  potential  defined  in 
equation  (2-1). 

i        +y\|  -1*  ,  the  imaginary  unit.  Also  used  as  an  index  in  the 
finite  difference  mesh  in  Section  IV. 

n        The  number  of  interior  points  in  the  finite  difference 
mesh  of  Section  IV.  Also  used  as  the  imaginary  part  of 
B,  i.e.  B=in. 

Re       Reynolds  number  based  on  the  mean  velocity  and  pipe  ra- 
dius. 


Time 


T,t,P,Q   Shorthand  notation  for  commonly  occurring  groups  of 
symbols  defined  in  equations  (2-18),  (2-19),  (3-U2), 
and  (3-U3). 

v        Velocity  vector  of  the  perturbation  flow  defined  by 
equation  (2-3). 

W        Complex  vector  potential  of  perturbation  velocity  defined 
by  equations  (2-1)  and  (2-2). 

x,r,9     Cylindrical  coordinates. 


A„  +  iAT  Complex  wave  number  of  the  perturbation  in  the 
x  direction. 


iBT  Complex  wave  number  of  the  perturbation  in  the  9  di- 
rection 

y  +  ^T  Complex  frequency  of  the  perturbation. 


r 


The  vorticity  transport  equation  (2-5)  expressed  in  ab- 
breviated notation  as  defined  by  equation  (2-6). 


P  '  r  >Vq  ^ie  components  of  p  in  cylindrical  coordinates  defined  by 
equation  (2-6). 

<3        Vorticity  vector  of  the  perturbation  flow  defined  in 
equation  (2-h). 


V  Linear  vector  operator  (nabla) . 

x  Vector  cross-product  operator. 

C  3  Brackets  enclosing  a  matrix. 

{  \  Brackets  enclosing  a  column  vector, 
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I.  INTRODUCTION 

The  research  reported  in  this  thesis  was  undertaken  in  an  attempt 
to  isolate  sources  of  error  causing  obviously  incorrect  computer  so- 
lutions for  the  three-dimensional  linearized  vorticity  transport  equa- 
tion for  pipe  Poiseuille  flow.  Solutions  obtained  using  the  theory  and 
program  developed  in  Ref.  1  had  indicated  decreasing  flow  stability 
with  decreasing  Reynolds  number,  a  result  which  is  clearly  inconsistent 
with  theory  and  experimental  data. 

Initial  analysis  of  that  program  confirmed  that  it  was  basically 
correct  as  presented  in  Ref.  1,  but  the  nature  of  the  soliztions  ob- 
tained indicated  one  or  more  errors  in  the  formulation  of  the  problem 
to  be  solved.   Analysis  of  this  problem  formulation  was  conducted  in 
three  phases.  The  first  phase  consisted  of  ensuring  accurate  expres- 
sion of  the  vorticity  transport  equation  in  terms  of  a  complex  velocity 
vector  potential.  The  second  and  most  significant  phase  of  the  analy- 
sis required  extensive  examination  of  the  conditions  required  to  satis- 
fy the  vorticity  transport  equation  at  the  singular  point  on  the  axis 
of  symmetry  of  the  pipe.  The  third  phase  entailed  establishing  a  fi- 
nite difference  approximation  of  this  equation  and  its  boundary  con- 
ditions for  adaptation  of  the  problem  to  a  standard  form  for  computer 
solution. 
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II.  THE  VORTICITY  TRANSPORT  EQUATION 

The  governing  equations  for  laminar  flow  of  an  incompressible  fluid 
with  constant  viscosity  are  the  continuity  equation  and  the  Navier- 
Stokes  equation.   Taking  the  curl  (vx)  of  the  latter  equation  and  intro- 
ducing a  perturbation  velocity  (v)  and  vorticity  (u>)  leads  to  the 
linearized  perturbation  vorticity  transport  equation  as  derived  in 
Appendix  A  of  Ref .  1 . 

Expressed  in  terms  of  the  complex  velocity  vector  potential,  W, 
which  has  the  form 

W(x,r,9,t)  =  (ixF(r)  +  erG(r)  +  iQH(r))eX  (2-1) 

where 

X  -  Ax  +  B9  ■»  Kt  (2-2) 

and 

v  =  vxW  (2-3) 

co  =  vxv  ,  (2-U) 

the  vorticity  transport  equation  becomes  three  simultaneous  fourth-order 


differential  equations  of  the  form 


DUF 


D3F 


2 

D  F 


l^\A  \  ♦  Cm31^d3g  I  +  <pg  +  *Dg)M 


|DUH 


3 


rDF^ 


(0^1  *»CN1"3)-i 


DG 
U>H 


- +  CDQ  +*DQ)iG 


(2-5) 


as  derived  in  Appendix  E  of  Ref.  1. 

Equations  (2-5)  may  be  expressed  in  the  abbreviated  form 
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n  V  =  <  o 


(2-6) 


0 


where  p  appears  to  be  a  set  of  three  coupled  equations  in  the  components 
of  W.  As  shown  in  Appendix  B  of  Ref.  1,  equations  (2-6)  actually  ex- 
press only  two  independent  conditions.  Therefore,  an  appropriate  linear 

combination  of  p  and  H  allows  expression  of  equations  (2-5)  as  a  set 
x      y 

of  two  equations  in  three  unknowns.  The  linear  combination 


±V    *   AH  =  0 


'x 


(2-7) 


does  not,  in  general,  uncouple  P  but  does  reduce  the  highest  order  deri- 
vative of  the  component  G(r)  in  equations  (2-5)  to  second  order. 

In  a  manner  similar  to  that  just  described  for  p,  Appendix  C  of 
Ref.  1  illustrates  the  redundancy  of  the  three  components  of  W.  This 
allows  arbitrary  selection  of  one  of  these  components  as  being  uniformly 
zero  for  all  r.  As  will  be  seen  later,  the  maximum  benefits  from  the 
consequences  of  the  linear  combination  expressed  by  equation  (2-7)  are 
obtained  if 

F(r)  =  0  .  (2-8) 

Incorporation  of  equations  (2-7)  and  (2-8)  in  equation  (2-5)  results 
in  expression  of  the  vorticity  transport  equation  in  the  form 

2r 


dCH..  ■  +i<j 


3  "\ 


D3H 


■  +  (D-1,1  -  *VW 


(cm;i-,d;d|j  ♦<£«;■] -»E»;i>-n 


(2-9) 


where 
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1X1 


Cm31  = 


0        -A 
Re 

0       0 

0        -2Al 

rRe 

0         B 

rRe 


(2-10) 


(2-11) 


oO  = 


dm 


-UAB 

Af  T 

r2Re 

V 

-t 

Re 

2B 
r2Re 

"ijAB 

3_  -  ■ 


(2-12) 


3 
r^Re 


(2-13) 


M 


(2-1U) 


PH"?})  t"  MM."  ?]"'}) 


[N21  .       0       -A 

jo       0 


"0 

^1 

r 

0 

B 

L_ 

r 

(2-15) 


(2-16) 


Tn'"]  =     T-2AB       A/J_ 

2  (      2 

r  V  r 


(7 "  j 


r 


B 


(2-17) 


II* 


and 


2    2 

t  2  A  +  B 


T  2  2A(1  -  r  )  -  t_ 

Re 


(2-18) 
(2-19) 


This  statement  of  the  vorticity  transport  equation  is  identical  to 
that  published  in  Ref.  1.  Note  that  equation  (2-8)  makes  it  unneces- 
sary to  include  F(r)  and  its  derivatives  in  equations  (2-9).  The  co- 
efficients of  this  function  and  its  derivatives  are  included  in  the 
arrays  represented  by  equations  (2-1 5)  through  (2-22)  of  Ref.  1  and 
have  been  verified.  The  reader's  attention  is  called  to  a  misprint 
ocurring  in  Ref.  1  concerning  the  development  to  this  point.  Equation 
(2-22)  of  Ref.  1  should  be  corrected  to  read  as  follows. 


m  •  r^t 


i 

i 


2AB 


'('"; 


^B 
2 


(2-20) 
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III.  BOUNDARY  AND  SPECIAL  CONDITIONS 

A.  BOUNDARY  CONDITIONS  AT  THE  WALL 

The  boundary  conditions  at  the  wall,  r=1,  are  derived  from  the  fact 
that  the  axial,  radial,  and  angular  components  of  the  flow  velocity  at 
the  wall  are  identically  zero.  These  components  are  obtained  respec- 
tively from  the  three  components  of  the  velocity  vector  resulting  from 
equation  (2-3). 

As  shown  in  part  two  of  Appendix  F  in  Ref  1  for  the  case  F(r)=0, 
these  components  imply  the  following  boundary  conditions  at  the  wall. 
G(1)  =  0  (3-la) 

H(1)  =  0  (3-1b) 

DH(1)  =  0  (3-1c) 

B.  CONDITIONS  ON  THE  AXIS  OF  SYMMETRY 

It  is  important  to  note  that  the  line  of  points  constituting  the  -  ~ 
axis  of  symmetry  is  not  a  boundary  but  rather  a  line  consisting  of  an 
infinite  number  of  singular  points  in  the  flow. 

It  should  be  noted  in  passing  that  all  functions  are  periodic  with 
respect  to  9.  Consequently,  as  shown  in  part  one  of  Appendix  G  in  Ref. 
1,  the  wave  number  of  the  perturbation,  B,  must  be  a  pure  imaginary 
number  such  that 

B  =  ni  ,    n  =  0,1,2,...  (3-2) 

Hence  references  to  B  in  this  paper  shall  be  understood  to  indicate  a 
pure  imaginary  number. 

The  primary  condition  to  be  imposed  is  that  the  vorticity  transport 
equation  be  satisfied  at  the  point  r=0.  Inspection  of  equations  (2-9) 
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through  (2-19)  shows  that,  in  general,  these  are  not  satisfied  for 
direct  substitution  of  r=0.  To  resolve  this  situation,  consider  the 
following  development. 

The  coefficients  of  the  functions  G(r)  and  H(r)  and  their  deriva- 
tives with  respect  to  r  up  to  fourth  order  in  equations  (2-9)  through 
(2-19)  may  be  expressed  in  the  following  form. 

oO-ovoi 

D>33  -LK3.j]l/r 


CM23  ■  CK2-2l1/r     *  £':2-ol  +  tM2+2l 


3 


C«ol  ■  [Vj,!1^  ♦£Mo-2l1/'2  +LK0-oV  UwV 

tN?l  ■  tN2-ol 

[No3  ■  [';o-2l1/r2  +  LHo-ol 


where 


[Vol 


M  " 


r- 

0 

-A 

Re 

0 

0 

0 


-2A 
Re 

B_ 
Re 


tM2-2l  - 


-hAB 
Re 

A  (3 

Re 

\je 

2B 

Re 

-   2B   ) 


*n 


(3-3) 
(3-U) 

(3-5) 
(3-6) 
(3-7) 
(3-8) 
(3-9) 
(3-10) 

(3-11) 
(3-12) 


(3-13) 
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CM2-ol  " 


[M2+2]  * 


IM1-31  = 


0 

2A2 

u 

"b 

-2A2" 

0 

_ 

0 

I  Re, 


0 


PS-0 


IxAB  A_(2B~   -   3; 

Re  Re 

if.  LJb2  - 1) 

Re  Re 


0  2A    f  1    -  A 


-A  -AB 

Re 


\         Re) 


CH1+11 


CMo-u) 


0       -2A 


tM0-2l 


0        2AB 
9 


-Uab(b 


Re 


"   +   1)       A_/3  +  B^(2  -  B^)> 
Re  I  ; 


-B2(B2   +1)  B    (B2   +    1) 


Ke 
UA2B 


Re 


tMo-o3  = 


l2AB^   ( 1    -  A_\     +  A_f       AB  /  3A  -   2\ 

\        Re;         Re  \Re  j 


(ir  -  1) 


(3-1U) 


(3-15) 


(3-16) 


(3-17) 


(3-18) 


(3-19) 


(3-20) 


(3-21) 


-Ua2b 


^(A2{2-R7V2B2) 

tMo+2~3  = 


A2/2(1    -  B2)    +  A2  (2 


-2AB 


0         -2A 


-2A^  0 


(3-22) 
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t",-il  ■ 
LNo-2"]  ■ 

[No-o]  ■ 


0 

-A 

0 

0 

"o 

-A" 

Q 

B 

P-2AB 

A(1 

-   B2) 

U 

B 

J 

r° 

3  n 

-AT 

U 

0 

(3-23) 


(3-2U) 


(3-25) 


(3-26) 


The  Maclaurin  series  representation  of  the  functions  G(r)  and  H(r) 
and  their  derivatives  may  be  expressed  according  to  the  following 
scheme. 


{ 
{ 


G" 

|G 

» 

■ 

H 

u 

J 

DG" 

- 

»      = 

DH 

2 
D  G 

■\           r 

-    — 

D2H 

J 

V 


G(r)  +  DG(r)  +  D2G(r)  +  D3G(r)  +  ...."" 

I(r)  +  DH(r)  +  D2H(r)  +  D3H(r)  +  

DG(r)  +  D2G(r)  +  D3G(r)  +  D^G(r)  + 

Drl(r)  +  D2H(r)  +  D3H(r)  +  D^H(r)  +  .... 

'D2G(r)  +  D3G(r)  +  B^Gtr)  +  D%(r)  +  . . . .' 

D2H(r)  +  D3H(r)  +  D^H(r)  +  D%l(r)  +  .... 


} 


where 


r 


{') 


^ 


2 
r  /2! 

r3/3! 

rVh! 


I 


(3-27) 


(3-28) 


(3-29) 


(3-30) 
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Substitution  of  equations  (3-3)  through  (3-30)  into  equation  (2-9) 
converts  the  vorticity  transport  equation  to  the  following  form. 


K-ol 


DaG(r)  +  ... 
DhH(r) 


D3G(r)  +  DUG(r)  +  ..7* 


1/r 


L  3*1^  D3H(r)  +  DUH( 


V. 


<Vr^ 


M 


D2G(r)  +  D3G(r)  +  D^G(r)  +  ..." 
D2H(r)  +  D3H(r)  +  D^H(r)  +  ... 


i 


1/r 
1/2 


^ 


tM2-ol 


D^G(r)  +  ... 
D2K(r)  +  ... 


-\ 


+  ["2  +  2"] 


D^G(r) 
>2T 


D  H(r)  +  ... 


tM1 -3l 


DG(r)  +  D2G(r)  +  D3G(r)  +  DijG(r)  +  ..."" 
DH(r)  +  D2H(r)  +  D3H(r)  +  D^H(r)  +  ... 


.DG(r)  +  D^G(r)  +  ... 
»DH(r)  +  D^H(r) 


"lX}(r)  + 
£M^lDH(r)  * 
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CMo-hl 


G(r)    +  DG(r)    +  D2G(r)    +  D3G(r)    +  D^r)    + 
_H(r)   +  DH(r)    +  D2H(r)   +  D3H(r)   +  D^H(r)    + 


r 


M 


G(r)    +  DG(r)    +  D^Gfr)   +   ... 
H(r)   +  DH(r)   +  D2H(r)   +   ... 


/r2^\ 


+  [%  q) 


-*[N2-o") 


G(r)   +   ..." 
H(r)    +   ... 

D2G(r)    +   ..7 
D2H(r)   +   ... 


.  \ +  M 


G(r)    +   . 
H(r)    +   ... 


•\  "^Ni-il  L»  + 


O/r 


'  *[N0-2] 


"   [No-ol 


G(r)    +  DG(r)    +  D^G(r)    +    ..." 
K(r)    +  DH(r)    +  D2H(r)    +    ... 


G(r)    +   ..." 
H(r)   +   ... 


< 


1/r 

1/2 


r 


v,  •     J 


1/r3 
1/2r2 
1/6r 
1/2U 


(3-31) 


Since  equations  (3-31)  must  be  satisfied  at  all  points  in  the  flow, 
they  must  be  satisfied  in  the  limit  as  r  approaches  zero.  A  method  of 
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determining  the  conditions  to  be  satisfied  at  the  singular  point,  r=0, 
is  described  below. 

Noting  that  as  r  approaches  zero  those  terms  in  equations  (3-31) 
containing  r  to  the  first  power  or  greater  may  be  neglected,  these 
equations  may  be  regrouped  as  a  power  series  in  r  and  expressed  in  the 
following  abbreviated  form. 

fo(or 

H(0), 

'G(O)' 
LH(0) 

'dg(o)T 


Cc,l 


r-h     ♦ 


D2G(0) 


tC2l 


*"  3'ilD2H( 


£Cbl 


-2 


D3G(0) 


/  fDhG(0)" 

V^Vh(o> 


^''mCojJ 

'D2G( 
LC8^D2H(0) 


-1 


where 


LSI  =  LVul 

•liABP 


Re 

2 

-B  P 

Re 


-AP(P  -  h) 

Re 

BP 
Re 


\ 


Cc2l  =  \.Mi-3l+tMo-al 

"-UB3       -AB2(P  -  $)' 


J 


Re 

Re 

-B1* 

2B 

Re 

Re 

'g(o) 

+  ^1h(0, 


(3-32) 


(3-33) 


(3-3U) 
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Cc3"3--[V2l+DIi-3l+(1/2)lMo-ul 


(3-35) 


-2ABP       -AF(P  -  k) 
Re  2Re 


2Re 


3BP 
2Re 


[ci3  ■  M  -tfcwl 


(3-36) 


2ABQ  AQ(P  -   2) 

2  2/2 

B  Q 


Ccsl 


+  ./T        -B/Q  -  a£\ 

Re  V  Re  / 

[vn+tM2-2i  +  (i/2tMi.3i  +  (i/6)tMo.j 


(3-37) 


f-2AB(P   +  3)        -AB2(Bli  +  U)" 


3Re 


6Re 


-B   (P  +  3)  2B_(P  +  3) 

Me  3Re 


[C6>[M1-ll    +[M0-2l    -*(Pl-ll    +LN0-21) 


(3-38) 


2ABQ       AB  Q 


[c7] 


B  Q         -2BQ 


rAB    (P  +  8)  -A(P  *  8)(P  +  h) 

^Re  2liRe 


-B*   (P  +  8) 
2hRe 


$P(P  +  8) 
2liRe 


J 


(3-39) 


IM 


LM2-oV  LMi-il +  (1/2)LMo-2l - 

*(lK2-ol   +  LH1-ll   +   (1/2'lK0-2l) 

AQ(P  +   2) 
2 
2' 


(3-1*0) 


B  Q  -   3A. 
Re 


■)    1(3Q  +  S 
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-   f-liAB2 


_ACAlQ  +  fl-2B2) 


(3-U1) 


(A  {Q  *£  }■ 2(p  -  2)j 


-2AE 


and 


P 

Q 


+   1 


2A   (\    -  k\    + 

V        He/ 


(3-1*2) 
(3-1*3) 


As  a  consequence  of  equations  (3-32),  the  following  conditions  must 
be  met  in  the  limit  as  r  approaches  zero. 


CciM 


G(0) 


VJ 


H(0) 


JDG(O) 

"■  2  1dh(o) 


D2G(0) 


1  'd'h(o) 


D3G(0) 


L  b J|  DJH(0) 


t^ 


DhG(0) 
J)4H(0) 


fD2G(0) 


G(0)' 
H(0) 


(3-1*1*) 


(3-1*5) 


(3-1*6) 


(3-1*7) 


(3-1*8) 


As  discussed  in  Section  II,  the  vorticity  transport  equation  as  ex- 
pressed by  equations  (2-9)  through  (2-19)  is,  in  general,  a  coupled  set 
of  differential  equations.  However,  if  B=0  these  equations  do  uncouple, 
For  this  special  case,  the  first  of  equations  (2-9)  becomes  a  homoge- 
neous fourth-order  differential  equation  in  H(r)  and  the  second  becomes 
a  homogeneous  second-order  differential  equation  in  G(r).  Thus,  the 


2U 


conditions  expressed  by  equations  ( 3— UU )  through  (3-U8)  may  now  be  ex- 
amined for  the  special  case  B=0  and  the  general  case  B>  0. 

1.  B=0 

Inspection  of  equations  (3-U5)  and  (3-U7)  shows  that  they  are 
identically  satisfied  for  all  values  of  the  applicable  odd-order 
derivatives  of  G(0)  and  H(0)  for  this  case.  Since  the  remaining 
equations  uncouple  for  this  case,  the  conditions  on  G(0)  and  H(0)  may 
be  studied  independently. 

a.  H(0) 

Sequential  inspection  of  the  first  of  equations  ( 3-l^U )  5 
(3-U6),  and  (3-I}8)  verifies  that  there  are  three  conditions  to  be  en- 
forced. These  three  conditions  are  as  follows. 

H(0)  =  0  (3-U9) 

D2H(0)  -  0  (3-50) 

ih(o)   =  o  (3-50 

b.  G(0) 

Although  the  second  of  equations  (3-Uh)  is  identically 
satisfied  for  arbitrary  G(0),  sequential  inspection  of  the  second  of 

equations  (3-U6)  and  ( 3-U8 )  yields  the  following  two  conditions  which 
must  be  enforced. 

G(0)  =  0  (3-52) 

D2G(0)  =  0  (3-53) 

2.  B>0 

With  the  exception  of  improbable  special  cases  for  which  the 
determinants  of  any  or  all  of  the  arrays  £C  "]  through  £C  "]  may  be  zero, 
the  conditions  to  be  met  for  this  case  are  as  follows. 
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G(0)  =  H(0)  =  0  (3-5U) 

DG(O)  =  DH(O)  ■  0  (3-55) 

D2G(0)  =  D2H(0)  =  0  (3-56) 

D3G(0)  =  D3H(0)  =  0  (3-57) 

D^G(O)  =  D^HCO)  =  0  (3-58) 

The  conditions  expressed  above  represent  a  marked  departure  from 
the  conditions  previously  thought  to  exist  at  r=0.  Before  the  situation 
was  properly  understood,  it  was  thought  that  the  required  conditions 
could  be  deduced  from  considerations  of  single-valuedness  at  this 
point.   This  approach  is  discussed  in  part  two  of  Appendix  G  in  Ref .  1 
but  must  now  be  discarded  as  insufficient.  In  the  interest  of  accuracy, 
the  reader's  attention  is  called  to  an  error  in  equation  (G-29)  of  that 
development.  That  equation  should  be  corrected  to  read  as  follows. 

v^O)  =  -iw^O)  (3-59) 
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IV.  NUMERICAL  METHODS 


Regrouping  equations  (2-9)  allows  expression  of  the  equations  to  be 


solved  in  the  following  format. 


c«;,i 


-nt»;i 


re 


lMol 


(U-1) 


The  coefficient  arrays  are  defined  by  equations  (2-10)  through  (2-19). 

As  discussed  in  Section  III,  the  nature  of  this  system  of  equations 
depends  on  the  value  of  the  complex  perturbation  wave  number,  B.  For 
the  special  case  B=0,  the  first  of  equations  (li-1)  becomes  a  homoge- 
neous fourth-order  differential  equation  in  H(r)  and  the  second  becomes 
a  homogeneous  second-order  differential  equation  in  G(r).  For  the 
general  case  of  B>0,  equations  (h-1)  do  not  uncouple  and  must  be  sol- 
ved simultaneously. 

An  additional  influence  of  B  is  reflected  in  the  character  of  the 
conditions  at  the  singular  point  r=0.   These  conditions,  as  developed  in 
Section  III,  are  summarized  below. 

For  the  special  case  B=0,  the  first  of  equations  (h-1)  must  simul- 
taneously satisfy  the  conditions 


H(0)  =  0 
2 


(U-2) 

(U-3) 
D4,H(0)  =  0  (U-U) 

and  the  second  of  equations  (h-1)  must  similarly  satisfy  the  following 
conditions. 


D  H(0)  =  0 
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G(0)  =  0  (U_5) 

D2G(0)  =  0  (U-6) 
For  the  general  case  where  B>0,  equations  (U-1)  must  simultaneously 
satisfy  the  conditions  belox^. 

G(0)  =  K(0)  =  0  (U-7) 

DG(O)  =  DH(O)  =  0  (^-8) 

D2G(0)  =  D2H(0)  =  0  (Ii-9) 

D3G(0)  =  D3K(0)  =  0  (U-10) 

dV.(o)  =  dSi(o)  =  0  •        (U-11) 

The  boundary  conditions  at  the  wall,  r=1,  are  a  consequence  only  of 
zero-velocity  viscous  effects  at  that  point  and  thus  do  not  vary  with  B. 
These  conditions  are  as  follows. 

G(1)  =  0  (U-12) 

H(1)  =  0  (U-13) 

DH(1)  -  0  (U-1h) 

Using  the  method  of  finite  differences,  the  functions  G(r)  and  H(r) 
may  be  approximated  by  a  finite  number  of  discrete,  evenly  spaced  un- 
knowns. As  shown  in  Figure  li-1  below,  the  non-dimensionalized  radius  of 
the  pipe  may  be  divided  into  a  one-dimensional  computational  mesh  of 
uniform  spacings  consisting  of  n  interior  points,  n+1  intervals,  and  n+2 
total  points,  including  the  boundary  point  at  r=1  and  the  singular 
point  at  r=0.  The  uniform  spacing  between  each  of  the  mesh  points  is 
%  as  defined  below. 

%  =  1/(n+1)  (U-15) 
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1 

r=o   - 

Figure  !i-1   Finite  Difference  Mesh 

By  taking  the  Taylor  series  expansion  of  the  function  H(r)  about 

the  i   mesh  point  in  terms  of  the  values  of  this  function  at  the  mesh 

points  i+2,  i+1,  i-1,  and  i-2  the  second-order  central  difference  ap- 

proximations  of  the  derivatives  of  H(r)  at  the  i   point  are  found  to 

have  the  following  form. 

DHi  =  (Hi+1  "  Hi-1)/2S  (li-l6) 

D2H.  =  (Hi+1  -  2H.  +  \^)/i2  (U-17) 

D3H.  -  (Hi+2  -  2(H1+1  -  H±-1)  -  H._2)/2S3  (U-18) 

D\  =  (H.+2  -  ltfli+1  +  6H.  -  lm±_,    +  H._2)/Sa  (11-19) 

where 

i  =  n, n-1, ...,3,2,1  (U-20) 
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The  approximations  for  the  derivatives  of  G(r)  are  determined  in  an 

identical  manner.  The  error  of  equations  (h-16)  through  ( U— 19)  is  of 

2 

the  order  of  magnitude  &  .   The  order  of  magnitude  of  this  error  may  be 

reduced  by  expanding  the  method  of  derivation  of  these  equations  to  in- 

J.  L 

elude  more  adjacent  points  on  either  side  of  the  i  ,  central,  point. 

It  is  important  to  note  a  discrepancy  appearing  in  Ref.  1  concerning 
the  development  to  this  point.  Careful  comparison  of  Figure  3-1  in  Ref. 
1  with  Figure  li-1  in  this  paper  reveals  that  the  direction  of  the 
labeling  schemes  used  to  depict  the  finite  difference  mesh  has  been  re- 
versed. Note  also  that  equations  (3-2)  in  Ref.  1  compare  exactly  with 
equations  (b-16)  through  (U- 19)  above.  V/hen  these  equations  are  used  in 
conjunction  with  the  labeling  scheme  used  in  Figure  3-1  in  Ref.  1,  the 
signs  of  the  odd-order  derivative . approximations  are  reversed.  This  was 
a  major  factor  in  producing  the  erroneous  solutions  obtained  from  the 
program  presented  in  that  paper. 

Substitution  of  the  central  difference  approximations  for  the  deri- 
vatives of  the  functions  G  and  H  in  each  of  equations  (Ii-1)  results  in 
a  set  of  n  linear,  algebraic  difference  equations  in  terms  of  the  un- 
known values  for  each  of  these  functions  at  each  of  the  n  interior 
points  of  the  finite  difference  mesh  depicted  in  Figure  li-1.  Since  each 
of  these  equations  is  of  the  form  of  a  linear  combination  of  the  i  , 
central,  mesh  point  and  the  two  adjacent  points  on  either  side  of  this 
central  point,  this  system  of  equations  consists  of  banded  coefficient 
arrays  multiplying  vectors  containing  the  unknown  values  of  the  func- 
tions at  each  of  the  n  interior  points.  By  using  this  technique  the 
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problem  is  converted  to  an  eigenvalue  problem  of  the  general  form 


tx){v}  -  *b1{v1  -  {0} 


U-21) 


with  the  basic  composition  of  the  arrays  \X]   and  QY}  and  the  vector  (vl 
as  illustrated  in  Figure  It- 2  below. 
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Figure  lt-2  Basic  Composition  of  the  Coefficient  Arrays  and 
the  Vector  of  Unknowns. 

The  exact  composition  of  £X~]   and  [_Y~\   depends  upon  the  value  of  B. 
These  arrays  are  established  by  the  subroutine  MSET2  in  conjunction 
with  the  function  subprograms  CHM1E1  and  CHM1E2,  which  compute  the 
numerical  value  for  each  element  in  these  arrays.  The  function  sub- 
program CHM1E1  provides  those  coefficients  required  from  the  first  of 
equations  (U-1)  and  CHM1E2  provides  those  coefficients  required  from 
the  second  of  these  equations.  For  the  special  case  B=0,  all  coeffi- 
cients contained  In  the  upper  right  and  lower  left  quadrants  of  Figure 
k-2   would  be  zero. 
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Of  particular  interest  are  the  difference  equations  whose  central 
points  are  adjacent  to  the  boundary  at  r=1,  point  n,    and  the  singular 
point  at  r=0,  point  1 .  To  amplify  this  matter,  consider  the  homoge- 
neous differential  equation  in  H(r)  resulting  from  setting  B=0  in 
equations  (U-1).  Use  of  equation. (U-1 9)  to  evaluate  'the  fourth  deri- 
vative of  this  function  at  point  n  of  the  finite  difference  mesh  re- 
veals that 

Similarly,  for  i=1  equation  (U-1 9)  becomes 

I)\    =  (H3  -  UH2  +  6H1  -  UHQ  +  H_1)Ah.  (U-23) 

Since  0-r-1,  equations  (U-22)  and  (U-23)  present  apparent  problems  in 

that  H  _  is  located  at  r=1  +  S  and  H  ,  is  located  at  r=-k,    i.e.  beyond 
n+2  -1 

the  allowable  range  of  values  for  r. 

These  inconsistencies  are  resolved  using  the  boundary  conditions 
expressed  by  equations  (U-1 2)  through  (U-1U)  and  the  conditions  at  r=0 
expressed  by  equations  (U-2)  through  (U-U).  Using  the  labeling  con- 
vention of  Figure  Lt —  1  and  equation  (U-1 3)  it  is  easily  confirmed  that 

H(1)  -  Hn+1  =  0  .  (U-2U) 

Comparison  of  equations  (U-1 U)  and  (U-1 6)  yields  the  following  relation- 
ships . 

DH(1)  =  DHn+1  -  (Hn+2  -  Hn)/2  =  0  (U-2*) 

H  +„  -  H  (U-26) 

n+2    n 

Equation  (U-26)  expresses  the  virtual  point  at  r=1  +  S  in  terms  of  the 
interior  point  n.  Equation  (U-22)  may  now  be  expressed  exclusively  in 
terms  of  boundary  and  interior  points  as  follows  (see  depiction  of 
equation  at  r=1-S  in  Figure  U-3  below). 
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Equation  at  r=1-S 
Equation  at  r=1-2S 
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Equation  at  r=1-5S 
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•   Central  point  of  the  difference  equation. 

0   Other  points  in  the  difference  equation. 

X   Boundary  or  singular  point  where  the  value  of  the 

variable  is  zero. 
D   Virtual  point  whose  coefficient  is  combined  with  the 

coefficient  of  the  point  indicated  by  an  arrow. 

Figure  U-3     Illustration  of  Basic  Method  Used  to  Include 
Virtual  Points  in  Difference  Equations. 


The  virtual  point  at  r=-S  can  be  expressed  in  terms  of  the  interior 
point  at  r=S  in  a  similar  manner.  Expressed  in  terms  of  the  labeling 
convention  of  Figure  U-1 ,    equation  (U-2)  becomes 

HQ  -  0  .  (U-27) 

Comparison  of  equations  (U-3)  and  (Ii-17)  shows  that 

(h-28) 
Direct  substitution  of  equation  (U— 27 )  in  equation  (U-28)  implies 

H_1  =  -H1  .  (U-29) 

This  leads  to  a  more  appropriate  expression  of  equation  (lj-23)  in  terms 


D2HQ  =   (H1    -   2HQ  +  H^/S2   =  0   . 


of  interior  points  as  follows. 

DItH1    =    (H3   -  UH?   +  6H1    -  liHQ  -  l^)/Sh 


(U-30) 


Consider  now  the  homogeneous  second-order  differential  equation  in 
G(r)  for  the  special  case  B=0.  Since  the  highest  order  derivative  of 
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this  function  contained  in  equations  (li-1 )  is  second-order,  inspection 
of  equations  (h-16)  and  (U-17)  verifies  the  fact  that  no  virtual  points 
are  required  in  the  vicinities  of  the  wall  or  the  axis  of  symmetry. 

As  previously  mentioned,  the  general  case  of  B>  0  requires  solution 
of  equations  (h-1)  as  a  system  of  coupled  equations.  Implementation  of 
the  conditions  at  r=0  for  this  problem  as  given  by  equations  (U-7) 
through  (U-11)  requires  special  attention. 

The  central  difference  approximation  of  the  fourth  derivative  of  H 
at  r=g  requires  a  virtual  point  at  r=-£.  Expression  of  this  virtual 
point  in  terms  of  the  interior  point  at  r=$  is  complicated  by  the  re- 
quirement that  the  five  conditions  for  the  function  H  given  by  equations 
(U-7)  through  (1|-11)  be  satisfied  simultaneously.  Using  equation  (U-16) 
to  satisfy  equation  (h-8)  seems  to  imply  that 

H_1  =  H1  .  (U-31) 

On  the  other  hand,  as  previously  shown  in  the  development  of  equation 

(U-29) 

H_1  =  -H.J  .  (U-32) 

This  contradiction  between  equations  (Ii-31 )  and  (U-32)  poses  a  pro- 
blem. However,  the  Maclaurin  series  approximation  for  the  function  H(r) 
and  its  derivatives  in  the  vicinity  of  r=0  is  helpful  in  resolving  this 
contradiction.  These  approximations  are  as  shown  below. 

H(r)  =  A  +  A.r  +  A  v_   +  A0r£  +  A,  r4  +  k^   *  0(r  ) 

2!  J?>\        FT   !TT* 


(U-33) 


DH(r)  =  A,  +  A  r  +  A  v_   +  A,  r£  +  A  r^  +  O(r^)  (h-3U) 

2!     3!     UT 


D2H(r)  =  A0  +  A_r  +  A,  r£  +  A^  +  0(r  ) 

1  ->  U9I       '-\\ 


3U 


(U-35) 


D3H(r)  =  A3  +  A^r  +  k^  *   0(r3)  (|t-36) 

D4H(r)  =  A^  +  A^r  +  0(r  )  (U-37) 

A  comparison  of  equations  (14-33)  through  (14-37)  with  equations  (Jb-7) 
through  (h-11)  indicates  that  at  r=0 

A0  =  A1  =  A2  =  A3  =  Aa  =  °  *  (l4_38) 

Substitution  of  equations  (a-38)  into  equation  (h-33)  yields  the  fol- 
lowing sixth-order  approximation  for  H(r)  in  the  limit  as  r  approaches 
zero. 

Evaluation  of  equation  (U-39)  for  r=+g  yields 


H(r)  -  k^P   +  0(rb)  (U-39) 


H(S)  -  A^  (U-UO) 

and 

h(-s)  =  a^x-hI  (U-Ul ) 

Inspection  of  equations  (h-hO)    and  (I4-I4I)  quickly  confirms  that  en- 
forcement of  equations  (It— 7 )  through  (li-11)  requires  that 

H_1  =  -H1  .  (1|-U2) 

Note  that  the  use  of  equation  (h-39)  as  an  approximation  for  the  func- 
tion H(r)  in  equation  (h-11)  does  not  alter  the  second-order  accuracy 
of  the  latter  approximation.  This  may  be  conceptualized  by  a  symbolic 
substitution  of  equation  (U-39)  in  equation  (U— 11)  as  shown  below. 

DhK1  =  ((H3  -  UH2  +  6H1  -  UHQ  +  H^)  +  0(*6))/gU  (U-U3) 

=  (H3  -  UH2  +  6H1  -  l4HQ  +  K_^)/lh   ♦  0(S2)  (U-UU) 
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With  the  formulation  of  the  problem  as  discussed  to  this  point 
accomplished,  i.e.  in  the  form  of  equation  (U-21),  the  method  used 
for  the  remainder  of  the  solution  is  summarized  in  the  following  steps: 
1 . )  Subroutine  CDMTIM  inverts  the  second  coefficient  array  in  equation 
(U-21),  array  £Y] .  CDMTIN  was  obtained  from  the  IBM  Library  routine 
CMTRIN  by  modifying  it  to  make  it  applicable  to  double  precision  arrays, 

2.)  Both  coefficient  arrays,  £X]  and  £Y]j  are  then  premultiplied  by 

-1 
£Y3   .   Since  multiplication  of  an  array  by  its  inverse  invariably  re- 

_i 
suits  in  the  identity  matrix,  £13  >   only  the  product  \y"\     £x3  is  com- 
puted using  subroutine  MULM.  This  converts  the  eigenvalue  problem  of 
equation  (ii-21 )  to  the  more  conventional  form 

(tzl  -itiDlv]  =  {o]  .  (a-l£) 

where 

tzl  =  W"1W  (U-Ii6) 

3.)  Since  all  the  programs  available  for  solving  equations  (h-h5)  re- 
quire that  the  real  and  imaginary  parts  of  the  elements  of  ^Z"]  be  pre- 
sented in  separate  arrays,  subroutine  DSPLIT  is  called  to  accomplish 
this. 

h»)      The  eigenvalues  of  equations  (1|-U£)  are  computed  by  subroutines 
EHESSC  and  ELRH1C  which  are  available  through  the  International  Mathe- 
matical and  Statistical  Library.  These  subroutines  reduce  the  matrix 
£Z"}  into  the  complex  upper  Hessenberg  form  and  then  solve  for  the 
eigenvalues.  The  results  are  then  passed  back  to  the  main  program,  R2. 
5.)  The  eigenvalues  thus  returned  are  listed  on  the  computer  output 
and  plotted  on  the  complex  plane  using  subroutine  PLOTP  from  the  IBM 
Library. 
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V.  RESULTS 

Owing  to  time  constraints ,  the  number  of  computer  solutions  obtained 
was  sufficient  only  to  confirm  the  general  validity  of  the  program  pre- 
sented in  this  paper.  For  this  reason,  a  detailed  discussion  of  flow 
stability  will  not  be  undertaken  here.  However,  the  following  brief 
discussion  should  provide  sufficient  information  for  interpreting  the 
data  presented  in  this  section. 

The  characteristics  of  the  flow  which  govern  the  stability  are  de- 
fined by  the  input  variables  A,  B,  and  Reynolds  number.  A  complete  set 
of  eigenvalues, 5f ,  is  obtained  for  each  chosen  combination  of  these 
variables.  Recalling  the  form  of  the  perturbation  velocity  vector  po- 
tential, VI,  as  presented  in  Section  II, 

W(x,r,0,t)  =  (e  F(r)  +  e  G(r)  +  e  H(r))exp(Ax  +  B9  +  *t)     (5-1) 
where 

it  is  readily  seen  that  positive  values  for  the  real  part,  tf~,  °f  the 
complex  perturbation  frequency  represent  an  exponential  growth  rate  in 
time  and  negative  values  of  Y  represent  an  exponential  decay  rate.  The 
algebraically  greatest  value  of  Y  contained  in  the  solution  set  is 
therefore  used  as  the  stability  criterion.  The  corresponding  eigenvalue 
is  said  to  be  the  least  stable  root.  The  flow  is  considered  to  be 
stable,  neutrally  stable,  or  unstable  with  respect  to  stationary  co- 
ordinates according  to  whether  the  least  stable  root  is  less  than,  equal 
to,  or  greater  than  zero,  respectively.  For  a  discussion  of  stability 
with  respect  to  moving  coordinates  see  part  D  of  Section  II  in  Ref.  1. 
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All  computations  vere  made  using  a  half  channel  mesh  containing  30 
interior  points  and  with  AT  held  constant  at  a  value  of  one.  Tables  I 
through  V  present  the  real  part  of  the  least  stable  root  for  each  of 
the  solutions  obtained  using  the  program  developed  in  this  paper.  In 
general,  these  solutions  are  encouraging  in  that  those  errors  which  had 
previously  caused  solutions  to  exhibit  decreasing  stability  with  de- 
creasing Reynolds  number  seem  to  have  been  resolved.  Results  appear  to 
be  reasonable  in  that  they  exhibit  stability  characteristics  which  are 
qualitatively  similar  to  those  obtained  for  the  plane  flow  problem  as 
detailed  in  Section  IV  of  Ref.  1.   Specifically,  for  the  range  of  Rey- 
nolds numbers  investigated,  all  flows  were  stable  for  Ao=0  while  insta- 

ft 

bilities  were  found  for  A=-0.05.  These  results  are  significant  in 
that  previous  investigations  of  fully  developed  pipe  flow  have  never 
accounted  for  the  instabilities  that  are  known  to  exist. 


Re 

V 

=  0 

AD  =  -0.05 
ft 

100 

-3.1*1* 

10"1 

-2.62  10"' 

200 

-2.88 

icf1 

-1.95  10"' 

5oo 

-1.80 

10-' 

-8.1i5  10~* 

1000 

-1.27 

10-' 

-2.99  10"* 

2000 

-8.93 

icf* 

8.57  10"^ 

5ooo 

-5.61 

10"  'd 

Ii.26  10"* 

10000 

-3.95 

icf* 

5.96  10"* 

Table  I   Y  D  f or  B,  =  0. 
ft     I 
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Re 

AR  =  0 

AR  =  -0.05 

100 

-3.05  10~1 

-2.31  10"1 

200 

-2.73  10*"1 

-2.02  10"1 

5oo 

-1.81*  10"1 

-9.31  10"2 

1000 

-1.33  10"' 

-3.98  10"2 

2000 

-9.58  10"2 

-3.91*  10~U 

5ooo 

-5.98  10"2 

3.7U  10~2 

1 0000 

-U.13  10"2 

5.68  10"2 

Table  II       tfR  for  B  =   1 


He 

AR  =  ° 

AR  =  -0.05 

100 

-h.91  10"1 

-3.38  10"1 

200 

-3.36  10"1 

-2.1*9  10"1 

5oo 

-2.09  10"1 

-1.19  10"1 

1000 

-1.1*8  10"' 

-5.1*8  10"2 

2000 

-1.05  io"1 

-9.5U  10"3 

5ooo 

-6.61*  10"2 

3.06  10~2 

10000 

-a. 71  io-2 

5.08  10'2 

Table  III       tfn  for  B  =   2 

K  1 
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Re 

ft 

■  0 

AR  -  -0.05 

100 

-U.91* 

10"1 

-1*.1*1*  io"1 

200 

-a.  27 

10"1 

-3.82  10"1 

500 

-2.89 

10"1 

-2.01  10"1 

1000 

-2.05 

10"1 

-1.13  10"1 

2000 

-1.U5 

io"1 

-5.05  io'2 

5ooo 

-9.16 

10"2 

1*.7U  10"J 

10000 

-6.1*9 

10"2 

3.25  10"2 

Table  IV      YD  for  B  =  3 
ft  x 


Re 

AR  = 

0 

AR  =  -0.05 

100 

-6.16 

10~1 

-5.70  io"1 

200 

-1*.69 

io"1 

-U..39  io"1 

500 

-3.1*1* 

10"1 

-2.87  10"1 

1000 

-2.75 

io"1 

-1.71*  10"1 

2000 

-1.87 

io"1 

-9.36  10~2 

5000 

-1.18 

10"' 

-2.25  10"2 

10000 

-8.36 

io"2 

1.33  10"2 

Table  V      V  D  f or  13  =  1* 
ft  x 
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